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Abstract 

A model-independent search for a narrow resonance produced in proton-proton colli¬ 
sions at ys = 8 TeV and decaying to a pair of 125 GeV Higgs bosons that in turn each 
decays into a bottom quark-antiquark pair is performed by the CMS experiment at 
the LHC. The analyzed data correspond to an integrated luminosity of 17.9 fb^^. No 
evidence for a signal is observed. Upper limits at a 95% confidence level on the pro¬ 
duction cross section for such a resonance, in the mass range from 270 to 1100 GeV, 
are reported. Using these results, a radion with decay constant of 1 TeV and mass 
from 300 to 1100 GeV, and a Kaluza-Klein graviton with mass from 380 to 830 GeV 
are excluded at a 95% confidence level. 
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1 Introduction 

Following the discovery of a Higgs boson (H) at the CERN LHC |[IH3l, with mass around 
125 GeV and properties so far consistent with the standard model (SM) of particle physics, it 
has become important to search for new resonances that decay into pairs of such Higgs bosons. 
While non-resonant pair production of the Higgs boson is allowed in the SM, the theoretical 
production cross section is approximately 10 fb |4J and well beyond the sensitivity of currently 
acquired data. However, several well-motivated hypotheses of physics beyond the standard 
model posit narrow-width resonances that decay into pairs of Higgs bosons, and could be pro¬ 
duced with large enough cross sections to be probed with existing data. The radion |5J and 
Kaluza-Klein (KK) gravitons in the Randall-Sundrum (RSI) |6| model of warped extra dimen¬ 
sions are examples of such resonances [7]. 

This letter reports the results of a model-independent search for the resonant pair production 
of Higgs bosons. The search for the narrow width resonance, denoted by X, is performed in the 
270-1100 GeV mass range. Data from proton-proton collisions at the LHG and recorded by the 
GMS experiment corresponding to an integrated luminosity of 17.9 ±0.5fb“^ at ^ = STeVis 
used. We perform this search for the case where both Higgs bosons decay into bottom quark- 
antiquark pairs (bb) [|8]. The main challenge of this search is to distinguish the signal of four 
bottom quarks in the final state that hadronize into jets (b jets) from the copious multijet back¬ 
ground described by quantum chromodynamics (QGD) in pp collisions. We address this chal¬ 
lenge by suitable event selection criteria that include dedicated b-jet identification techniques 
and a model of the multijet background that is validated in data control regions. Our results 
may be compared with a search performed by the ATLAS experiment [93 that also probes the 
physics of resonant Higgs boson pair production, albeit in the channel where one Higgs boson 
decays to bottom quarks and the other decays to photons. 

2 Detector and Event Reconstruction 

A detailed description of the GMS detector, together with a description of the coordinate system 
used and the relevant kinematic variables, can be found in Ref. IllOll . The central feature of the 
GMS apparatus is a superconducting solenoid of 6 m internal diameter that generates an axial 
magnetic field of 3.8 T. Within the solenoid volume are a silicon pixel and strip tracker, a lead 
tungstate crystal electromagnetic calorimeter, and a brass and scintillator hadron calorimeter. 
Muons are detected and their properties measured in gas-ionization detectors embedded in a 
steel flux-return yoke outside the solenoid. Jets are reconstructed using the anti-kj clustering 
algorithm IHTIIT^ with a distance parameter of 0.5 applied on the collection of particle candi¬ 
dates reconstructed by the particle-flow (PF) algorithm The PF algorithm reconstructs 

and identifies each individual particle with a combination of information from the various el¬ 
ements of the GMS detector. To mitigate the effect of additional particles that do not originate 
from the hard interaction in jet reconstruction, we subtract charged hadrons that do not arise 
from the primary vertex associated with the jet from the collection of clustered particles. Fur¬ 
ther, an average neutral energy density from particles not arising from the primary vertex is 
evaluated and subtracted from the jets fl5l . Energy corrections for the jets are determined as 
functions of the jet transverse momentum pj and pseudorapidity |//|. Jet identification cri¬ 
teria Kl^ to reject detector noise misidentified as jets, and jets not originating from the hard 
interaction are also applied. 

In order to identify (tag) b jets, we rely on the fact that bottom quarks hadronize into b hadrons 
which have decay lengths of the order of cr = 450 pm. Thus, their decay products originate 
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4 Event Selection 


from secondary vertices made of tracks that have impact parameters with respect to the pri¬ 
mary vertex of a similar scale. The pixel tracker provides an impact parameter resolution of 
about 15f/m for charged tracks with |)y| < 2.4. To maximize the b-tagging performance of 
the detector, we combine the output discriminants of several b-tagging algorithms described 
in Ref. flTl with a trained artificial neural network. This we call the combined multivariate 
(CMVA) algorithm. In particular, we combine the outputs of the combined secondary vertex 
(CSV) tagger that uses secondary vertices identified by the inclusive vertex finder (IVF) algo¬ 
rithm HSl, the jet probability (JP) tagger, and the two soft lepton taggers. 

The first level of the trigger, consisting of customized processors, collects data for this analysis 
using information from the calorimeters and requires two jets to exceed pj thresholds of 56 
or 64GeV, depending on luminosity conditions. The second level of the trigger, consisting 
of software algorithms executed on a farm of commercial processors, uses information from 
the entire detector to reconstruct PF jets, and requires four PF jets with \r]\ < 2.4 and pj > 
30 GeV, of which two jets must have pj > 80 GeV. Further, to record signal events and reject 
background QGD multijet events, two jets are required to be tagged by the GSV b-tagging 
algorithm implemented at the trigger. 


3 Simulated Samples 

To model the production of a generic narrow-width spin-0 resonance, we use a Monte Garlo 
simulation of the RSI radion produced through gluon fusion. The angular distributions of a 
spin-2 resonance are distinct from those of a spin-0 resonance, and result in different kinematic 
distributions. Therefore, we evaluate the signal efficiencies for a narrow-width spin-2 reso¬ 
nance from a separate simulation of the first excitation of the KK graviton produced through 
gluon fusion in the same extra dimension scenario as the radion. The resonance is forced to 
decay to a pair of Higgs bosons where both Higgs bosons decay to bb. Samples of these signal 
events, as well as background events from diboson, W+jets, Z+jets and top-quark pair pro¬ 
duction (tt) processes, are generated using the MadGraph 5.1 |fl9 ll program interfaced with 
PYTHIA 6.4 ||20| for parton showering and hadronization. QGD multijet event samples are sim¬ 
ulated with the PYTHIA 6.4 program. A sample of events where the Higgs boson is produced in 
association with a Z boson is simulated using the powheg event generator ||2TI - |2^ interfaced 
with the HERWIG++ 1241 program for showering and hadronization. We set the PYTHIA 6.4 
parameters for the underlying event to the Z2* tune [23. The response of the GMS detector is 
modeled using Geant4 [261 . 

On average, 21 pp interactions occurred per bunch crossing in the data used in this analysis. 
Additional simulated pp interactions overlapping with the event of interest were added to the 
simulated samples to reproduce the distribution of the number of primary vertices per event 
reconstructed in data. 

4 Event Selection 

The trigger-level jet pj thresholds confine our search for a narrow-width X —I HH —I bbbb 
resonance to masses above 270 GeV. Beyond nix ~ 800 GeV, the selection efficiency is increas¬ 
ingly limited by the merging of jets from the same Higgs boson, and we curtail this search at 
1100 GeV. The kinematic distributions of the decay products vary substantially over this mass 
range. Therefore, to optimize the search sensitivity, we use different event selection criteria 
in three main kinematic regions: the low-mass region (LMR) for mass hypotheses from 270 to 
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Figure 1: Illustration of the SR, SB, VR and VS kinematic regions in the nin^) plane used to 

motivate and validate the parametric model for the QCD multijet background. The quantities 
?«Hi and whz are the two reconstructed Higgs boson masses. The distribution in data events 
after b-tagging and kinematic selections is shown with the SR blinded. 

450 GeV, the medium-mass region (MMR) for masses from 450 to 730 GeV, and the high-mass 
region (HMR) for masses from 730 to 1100 GeV. 

Event selection begins with the identification of events containing at least four jets in the cen¬ 
tral region of the detector (|?/| < 2.4) that are b-tagged and have px > 40 GeV. To b-tag a jet, 
we require it to pass a working point for the GMVA algorithm that maximizes the sensitiv¬ 
ity of this search. For jets with px > 40 GeV and | p | < 2.4 this working point yields a 75% 
efficiency for tagging jets originating from b hadrons and a mistagging rate of 3% for light- 
flavor jets. For the LMR, we combine these b jets into pairs to create HH candidates such that 
|?mh ~ 125 GeV I < 35 GeV for each candidate Higgs boson. The mass resolution on the Higgs 
boson in the LMR is found to be approximately 9 GeV. Selected HH candidates are required to 
have at least two jets with px > 90 GeV. In the MMR, signal events have large Lorentz factors 
for the Higgs boson candidates. Therefore, HH candidates for this region are constructed from 
four jets such that the AR = VAp^ -|- A(^^ between the jets associated with an H candidate re¬ 
main within 1.5, where Ap and A^ are the differences in the pseudorapidities and azimuthal 
angles of the two jets. For the HMR, we use the same criteria used in the MMR with an ad¬ 
ditional requirement of px > 300 GeV on one of the H candidates to better discriminate signal 
events from background. In all three regions, in case of multiple HH candidates in an event, the 
combination with the smallest |?«Hi ~ ?^H 2 l is chosen. Having identified the two Higgs boson 
candidates in each event, we plot their masses, mHj and niYi^, on a two-dimensional histogram 
as shown in Fig. Hi and H 2 are chosen at random from the two reconstructed H candidates. 
As the final selection criterion applied in each of the three mass hypothesis regions, we re¬ 
quire events to fall within the signal region (SR) defined as VAm^^ -|- Am^^ < 17.5 GeV, where 
A?nHi 2 = ?”Hi 2 “ 125 GeV. 
































4 


5 Signal Modeling 


Table 1: Efficiencies of the event selection criteria for generic spin-0 and spin-2 resonances 
decaying to a pair of Higgs bosons in the four b jet final state at representative masses. 


selection eff. (%) "^oo* 

spin-0 0.05 

spin-2 0.09 


Mass (GeV) 
W'TOO900”TT00 

2.3 4.9 4.6 2.2 

3.4 6.6 5.4 2.3 


(8 TeV) 



Figure 2: A maximum likelihood fit to the nix distribution of simulated signal events for the 
700 GeV mass hypothesis. The distribution is fitted to a Gaussian core smoothly extended on 
both sides to exponential tails. Here n is the number of degrees of freedom in the fit. 

The efficiencies of these selection criteria for spin-0 and spin-2 resonances at representative 
masses are shown in Table The major loss in efficiency for all mass hypotheses comes from 
the b-tagging requirement for 4 jets. For the 300 GeV mass hypothesis, this is compounded by 
the trigger inefficiency. The distribution of the aforementioned AR between jets from a single 
Higgs boson is narrower for the spin-2 resonance, and thus requiring AR < 1.5 results in a 
higher efficiency for it. 


5 Signal Modeling 

For signal events, the aforementioned event selection criteria are expected to produce a sharp 
peak in the nix distribution over a relatively featureless background from events arising from 
SM processes. The interference between SM background processes and the narrow resonant 
signal is expected to be negligible. To search for signal events at various mass hypotheses, we 
fit the mx distribution in data events in the SR to a parametric model for the signal peak on 
top of parametric models appropriate for components of the SM background. This procedure 
is performed for the FMR, the MMR, and the HMR separately. 

To improve the mass resolution of the signal X —> HH resonance, we perform a fit that con¬ 
strains the invariant masses of the Higgs boson candidates. In the fit, the momenta of the re- 
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constructed b jets are allowed to float within their expected resolutions. Since the uncertainty 
in the reconstructed mass of the Higgs boson candidate due to the measurement of jet direc¬ 
tion is smaller than that due to the measurement of jet energy, this constraint mainly affects 
the latter. This fit improves the invariant mass resolution of the reconstructed signal resonance 
by 20^0%, depending on the mass hypothesis. Extensive tests in background-dominated con¬ 
trol regions in data show that no artificial structures are introduced in the background mass 
distributions by this procedure. 

We build the parametric model for each signal mass hypothesis by fitting the shape of the mx 
distribution of simulated events that are accepted by the selection criteria and corrected for dif¬ 
ferences between data and simulation. A sum of two Gaussian functions, requiring five param¬ 
eters, is used for the LMR fit to account for tails in the distribution from incorrect combinations 
of jets. In the MMR and the HMR, we fit a function with a Gaussian core smoothly extended 
on both sides to exponential tails, such that the function is continuous both in its value and its 
first derivative. This requires two parameters for the mean and width of the Gaussian function, 
and two other parameters for the exponential tails on both sides. An example of a paramet¬ 
ric model for the MMR signal obtained through this procedure is shown in Fig. While the 
model is constructed for the mass hypothesis of 700 GeV, its Gaussian core peaks at 714 GeV 
and has a width of 21 GeV. This mass shift is found to be linear in nix and occurs due to the 
aforementioned constraint of jet momenta to Wh- 


6 Background Modeling 

While the composition of background events in the SR is expected to be dominated by QGD 
multijet processes, we find through simulation that tt production contributes approximately 
22%, 27%, and 24% in the LMR, MMR, and HMR, respectively. We also find that Z+jets, ZZ, 
and ZH processes contribute less than 1% of the background and therefore neglect them in this 
analysis. The nix distribution of these tf events is found to be somewhat different in shape from 
that of QGD multijet events, and therefore we treat it as a distinct component of the background 
and model it with a parametric form. We obtain this parametric form by fitting the shape of 
the mx distribution of simulated tf events accepted by the event selection criteria to a function 
with a Gaussian core smoothly extended to an exponential tail on the high side. This function, 
henceforth referred to as GaussExp, is continuous in its value and its first derivative. It has two 
parameters for the mean and width of the Gaussian function and one parameter for the decay 
constant of the exponential tail. This model is normalized to a tf cross section of 234 pb |[27l , and 
is allowed to float with a systematic uncertainty of 15% in the final fit to account for theoretical 
and measurement uncertainties in our kinematically boosted region. 

We use the GaussExp parametric model to fit the mx distribution of the QGD multijet com¬ 
ponent of the background in the SR. With the SR kept blinded, we motivate and validate this 
choice of parametric model by the fact that it fits well the shape of the mx QGD multijet back¬ 
ground distributions in several different regions of the {mn ^, ma ^) plane depicted in Fig. [^and 
described below. We do not aim to predict the parameters of the model in the SR from the other 
regions. These fits are performed for the LMR between 260 and 650 GeV, for the MMR between 
400 and 900 GeV, and for the HMR between 600 and 1200 GeV. In each case the tf contribution, 
as expected from simulation, is subtracted. 

We define a sideband region (SB) to the SR as 17.5 GeV < V -|- < 35 GeV and 

Aj^Hj < 0. For events in this region, the mx distribution is expected to be kinemati¬ 
cally similar to that for events in the SR, since in each of the sidebands one of the reconstructed 
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6 Background Modeling 



17.9 fb ' (8 TeV) 





Figure 3: The mx distributions in data (after the tf background has been subtracted) in the SB 
of the MMR (top left), the CR of the MMR (top right), the VS of the LMR (bottom left), and the 
VR of the LMR (bottom right). The distributions are fitted to the GaussExp parametric model. 
The shaded regions correspond to ±1(7 variations of this fit. Here n is the number of degrees 
of freedom in each fit. The pull, for a given bin, is defined as the number of data events minus 
the value of the fit model, divided by the uncertainty in the number of data events. 
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Higgs boson masses is slightly higher in value than for events in the SR while the other is 
slightly lower. As an example, Fig.|^top left shows the fit performed for events in the SB pass¬ 
ing the HMR selection. Another set of events that pass the kinematic requirements of the event 
selection criteria in the SR region of the (wjhj/ ^Hz) plane but required to have one of the four 
jets not be b-tagged is selected to further test the applicability of the GaussExp model in de¬ 
scribing the nix distribution of the QCD multijet background in a different but kinematically 
similar region. This is called the data Control Region (CR), and the fit for these events, that 
would have otherwise passed the HMR selection, is also shown in Fig. on the top-right. In 
both cases, the goodness of the fit, characterized by the per degree of freedom, is found to 
be reasonable. 

These two cases already lend significant confidence to the choice of the GaussExp parametric 
model for the SR. However, we carry out further checks in neighboring validation regions (VR) 
with a corresponding sideband (VS) that are defined similarly to the SR and SB regions but with 
njHi = wJHz centered at different values. The good fits for the nix distributions in these regions 
not only demonstrate the applicability of the GaussExp model to describe these kinematically 
distinct QCD multijet events, but also that events in the VR are in fact kinematically similar 
to those in the VS. As examples, Eig. |^bottom-left and bottom-right plots show the results of 
these fits for the LMR selection for the VS and the VR, respectively, both centered at Whi = 
= 90 GeV. We obtain similar results for the VR centered at Whi = w^Hz = 107.5,142.5 and 
160 GeV. 

While the GaussExp function fits well the nix distribution from QCD multijet events in all these 
distinct regions and therefore can be expected to be a good approximation of the parametric 
form of the true parent distribution for events in the SR, other similar parametric models could 
be chosen instead. Therefore, a systematic uncertainty associated with the choice of this para¬ 
metric model is evaluated by assuming a 7th order Bernstein polynomial, which also fits the 
nix distribution well in the SB, to be the true distribution. Pseudo-datasets are generated from 
this polynomial function and fitted with the GaussExp function as well as other polynomial 
functions to compute biases in the reconstructed signal strength. This procedure is performed 
for each mass hypothesis. These biases are found to be of the order of 100 fb for the LMR, 
10 fb for the MMR, and 20 fb for the HMR. We account for this bias as a signal-shaped system¬ 
atic uncertainty in the background model with normalization centered at zero and a Gaussian 
uncertainty with standard deviation equal to the bias. 


Table 2: Relative systematic uncertainties in the selection efficiencies for signal and tt events in 
the LMR, the MMR, and the HMR. 


Source of systematic 
uncertainty 

Impact in 

LMR (%) 

Impact in 

MMR (%) 

Impact in 

HMR (%) 

Signal 

tf 

Signal 

tf 

Signal 

tf 

Jet energy scale 

0.1-0.2 

0.8 

0.0-0.2 

0.1 

0.1-0.3 

1.9 

Jet energy resolution 

2.4-7.0 

2.7 

5.5-7.0 

2.2 

4.9-5.3 

7.2 

b-tagging scale factor 

12.7 

12.7 

12.7 

12.7 

12.7 

12.7 

Trigger scale factor 

6.0-18.8 

9.1 

6.1-8.0 

7.2 

7.4-7.9 

7.2 


7 Systematic Uncertainties 

Sources of systematic uncertainties that affect the selection efficiencies for signal and tf events 
are listed in Table |^and described below. We vary the jet energy scale |28 | by its uncertainty as 
a function of jet px and tj and find that this affects signal efficiencies by a relative factor of up 
to 0.2% and tf efficiencies by up to 0.8%. We evaluate the effect of uncertainty in the jet energy 
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8 Results 


resolution by varying the jet energies according to the measured uncertainty This is found to 
affect signal efficiencies by 2-7%, and tf efficiencies by 2-3%. These uncertainties affect not 
only the normalizations but also the parameters of the signal and tf models, and are taken into 
account as nuisance parameters in the final fit. 

The trigger efficiencies for signal and tf events are evaluated approximately by passing gen¬ 
erated events through a trigger simulation. We then correct these efficiencies for differences 
between simulation and data by computing the difference in a tf-enriched control region ob¬ 
tained using a trigger that requires at least one muon with pj > 24 GeV. Further, event selection 
criteria requiring at least one muon with pj > 40 GeV and at least four jets in the central region 
of the detector with pj > 40 GeV are applied. The data-to-simulation correction factor is char¬ 
acterized by the pj and GMVA discriminants of the relevant jets. Uncertainties in this factor 
impact signal efficiencies by 6-18%, and tf efficiencies by 7-9%. 

The b-tagging efficiencies of the GMVA algorithm for signal and tf events are also evaluated 
approximately through simulation and then corrected by a data-to-simulation comparison. The 
comparison is performed in the same tf-enriched control region as the calculation of the trigger 
efficiency. The correction factor for the b-tagging efficiencies is consistent with unity. The 
uncertainty in this factor for four b jets is evaluated to be 12.7%. 

Additionally, the yields of signal events for a given production cross section and tf events are 
both affected by a 2.6% uncertainty in the measurement of the integrated luminosity 129 1. 


8 Results 

The nix distribution that we observe in data within the SR, along with a binned maximum- 
likelihood fit with the aforementioned parametric background models, are shown in Fig. We 
compute the observed and expected upper limits on the cross section for pp —> X —> HH —> 
bbbb at a 95% confidence level (GL) using the modified frequentist GLg method Il30l l3Tl by 
fitting the data with the parametric signal, tf, and QGD multijet models. This is done separately 
in the disjoint ranges of nix for the individual regions described in Section]^ and the limits are 
presented together in Fig. These limits are shown for the spin-0 resonance on the left, and the 
spin-2 resonance on the right. The green (dark) and yellow (light) bands respectively represent 
the Ic and 2ct confidence intervals around the expected limits. The observed upper limits 
lie within 2(7 of the expected upper limits, and thus we conclude that there is no significant 
deviation from the background-only hypothesis. 

The theoretical cross section for the production via gluon fusion of a radion that decays to a pair 
of Higgs bosons |[32l that each in turn decays to a bb pair with a branching fraction of 58% 1331 
is calculated using MadGraph 5.1 Il34l and superimposed on the experimental cross section 
limit for the spin-0 resonance in the plot on the left. In this calculation, the correction factor 
used to account for next-to-leading-order effects for electroweak couplings Il35l and next-to- 
next-to-leading-order effects for QGD couplings l36l is identical to that used for Higgs boson 
production through gluon fusion. The warped extra dimension scenario for this radion has the 
product of the curvature, k, and half the circumference of the extra dimension, L, set to 35, a 
radion decay constant of Ar = 1 TeV, and no radion-Higgs boson mixing. The theoretical cross 
section for the radion has an uncertainty of approximately 15% that is not used to compute 
the experimental limits on spin-0 resonance production shown in Fig.|^ Masses for the radion 
between 300 and 1100 GeV are excluded at a 95% GL. A similarly calculated theoretical cross 
section for the KK graviton as the resonance X, in the same warped extra dimension scenario, 
is overlaid on the limit for the spin-2 resonance in the plot on the right. Masses for such a 
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Figure 4: The nix distribution in data in the SR between 260 and 650 GeV of the LMR (top left), 
between 400 and 900 GeV of the MMR (top right), and between 600 and 1200 GeV in the HMR 
(bottom). All distributions are fitted to the background-only hypothesis for illustration, show¬ 
ing the relative contributions of the QGD multijet (dashed-dotted red) and tt (dashed green) 
processes. The pull, for a given bin, is defined as the number of data events minus the value of 
the background-only fit, divided by the uncertainty in the number of data events. Also for il¬ 
lustration, we overlay the signal models of the spin-0 resonance (dotted blue) corresponding to 
mass hypotheses and production cross sections of 350 GeV and 653 fb for the LMR, 700 GeV and 
17.6 fb for the MMR, and 900 GeV and 8.1 fb for the HMR. These cross sections correspond to 
the observed upper limits, which are computed for signal mass hypotheses from 270 to 450 GeV 
in the LMR, from 450 to 730 GeV in the MMR, and from 730 to 1100 GeV in the HMR. 
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9 Summary 
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Figure 5: The observed and expected upper limits on the cross section for pp —> X ^ HH —> 
bbbb at a 95% confidence level, where the resonance X has spin-0 (left) and spin-2 (right). 
The theoretical cross section for the RSI radion, with Ar =1 TeV, kL = 35, and no radion- 
Higgs boson mixing, decaying to four b jets via Higgs bosons is overlaid on the left plot. The 
theoretical cross section for the first excitation of the KK-graviton for the same parameters is 
overlaid on the right plot. 


graviton are excluded at a 95% CL between 380 and 830 GeV. 


9 Summary 

We have presented a model-independent search by the CMS experiment at the LHC for a nar¬ 
row resonance produced in proton-proton collisions at \/s = 8 TeV and decaying to a pair of 
125 GeV Higgs bosons that in turn each decays into a bottom quark-antiquark pair. The an¬ 
alyzed data correspond to an integrated luminosity of 17.9 fb^^. No evidence for a signal is 
observed. Upper limits at a 95% CL on the production cross section for such sprn-O and spin-2 
resonances, in the mass range from 270 to 1100 GeV, are reported. Using these results, a radion 
with decay constant of 1 TeV and mass from 300 to 1100 GeV, and a Kaluza-Klein graviton with 
mass from 380 to 830 GeV are excluded at a 95% confidence level. 
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